knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, fig.align='center')
source("./packages.R")
source("./read_data.R")
Mothers & Newborns
- 5 phenols
- 3 parabens
- 9 phthalates
Outliers
mod = lm(WISC ~ P1 + SEX*P1 + SEX + M_IQ + ALCOHOL + M_EDU + M_AGE +
MARITAL_STATUS + HOME_SCORE + mat_hard
, data = whole)
summary(mod)
##
## Call:
## lm(formula = WISC ~ P1 + SEX * P1 + SEX + M_IQ + ALCOHOL + M_EDU +
## M_AGE + MARITAL_STATUS + HOME_SCORE + mat_hard, data = whole)
##
## Residuals:
## Min 1Q Median 3Q Max
## -43.649 -8.209 0.684 8.871 28.071
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 59.92308 9.25283 6.476 0.000000000425 ***
## P1 -0.32796 1.04128 -0.315 0.753029
## SEXFemale -1.79678 2.28932 -0.785 0.433208
## M_IQ 0.22730 0.06017 3.777 0.000194 ***
## ALCOHOLYes -5.78965 1.70620 -3.393 0.000791 ***
## M_EDU< High school degree -1.10596 1.65639 -0.668 0.504883
## M_AGE 0.24590 0.17339 1.418 0.157266
## MARITAL_STATUSEver Married -0.02059 1.70476 -0.012 0.990370
## HOME_SCORE 0.40270 0.12361 3.258 0.001262 **
## mat_hardYes -2.23735 1.53031 -1.462 0.144864
## P1:SEXFemale -0.15397 1.40791 -0.109 0.912998
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 12.45 on 278 degrees of freedom
## (54 observations deleted due to missingness)
## Multiple R-squared: 0.164, Adjusted R-squared: 0.134
## F-statistic: 5.455 on 10 and 278 DF, p-value: 0.0000002148
plot(mod)




# 274 influential based on Cooks distance, residual vs. leverage plot
ols_plot_cooksd_chart(mod)

ols_plot_resid_lev(mod)

# get rid of highest
whole_m1 = whole[-which.max(whole$P1),]
mod2 = lm(WISC ~ P1 + SEX*P1 + SEX + SEX + M_IQ + ALCOHOL + M_EDU + M_AGE +
MARITAL_STATUS + HOME_SCORE + mat_hard
, data = whole_m1)
summary(mod2)
##
## Call:
## lm(formula = WISC ~ P1 + SEX * P1 + SEX + SEX + M_IQ + ALCOHOL +
## M_EDU + M_AGE + MARITAL_STATUS + HOME_SCORE + mat_hard, data = whole_m1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -43.703 -8.060 0.495 9.058 28.178
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 61.28437 9.24589 6.628 0.000000000177 ***
## P1 -0.36972 1.03731 -0.356 0.721795
## SEXFemale -0.25220 2.43462 -0.104 0.917572
## M_IQ 0.22489 0.05994 3.752 0.000214 ***
## ALCOHOLYes -5.62428 1.70172 -3.305 0.001075 **
## M_EDU< High school degree -1.39491 1.65737 -0.842 0.400716
## M_AGE 0.24692 0.17269 1.430 0.153885
## MARITAL_STATUSEver Married -0.33112 1.70649 -0.194 0.846292
## HOME_SCORE 0.38196 0.12364 3.089 0.002210 **
## mat_hardYes -2.57323 1.53536 -1.676 0.094870 .
## P1:SEXFemale -1.61214 1.61733 -0.997 0.319735
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 12.4 on 277 degrees of freedom
## (54 observations deleted due to missingness)
## Multiple R-squared: 0.1698, Adjusted R-squared: 0.1398
## F-statistic: 5.666 on 10 and 277 DF, p-value: 0.0000001013
plot(mod2)




# cooks dist for 264
ols_plot_cooksd_chart(mod2)

ols_plot_resid_lev(mod2)

whole[c(319, 307),]
| SID | ETH | M_EDU | MARITAL_STATUS | SMOKER_IN_HOME | SEX | HOME_SCORE | ALCOHOL | M_AGE | M_IQ | GEST | mat_hard | WISC | P1 | P2 | WISC_VC | WISC_PR | WISC_PS | WISC_WM |
| 1229 | African American | < High school degree | Ever Married | No | Female | 40 | No | 30 | 102 | 40 | Yes | 113 | 10.1 | 1.99 | 93 | 125 | 115 | 110 |
| 1209 | African American | ≥ High school degree or equivalent | Never Married | Yes | Female | 36 | No | 36.1 | 80 | 39 | No | 103 | 7.14 | 2.57 | 106 | 106 | 83 | 107 |
max(whole$P1) # 319
## [1] 10.0907
grubbs.test(log(whole$P1)) # highest values is an outlier
##
## Grubbs test for one outlier
##
## data: log(whole$P1)
## G = 3.88600, U = 0.95572, p-value = 0.01474
## alternative hypothesis: highest value 2.31161455316489 is an outlier
grubbs.test(log(whole_m1$P1))
##
## Grubbs test for one outlier
##
## data: log(whole_m1$P1)
## G = 3.3899, U = 0.9662, p-value = 0.1084
## alternative hypothesis: highest value 1.96633990278861 is an outlier
Demographics
# Ran BN2MF on this subset
ppp_343 = left_join(ppp_cov, mn_outcome) %>% mutate(Subsample = "BN2MF Subsample") %>% rename(HOME_SCORE = HOMETOT)
# this subset had WISC measured
ppp_311 = ppp_343 %>% filter(!is.na(WISC)) %>% mutate(Subsample = "WISC Subsample")
# this subset had all covariates measured
ppp_289 = ppp_311 %>% filter(!is.na(M_IQ) & !is.na(HOME_SCORE)) %>% mutate(Subsample = "Complete Cases")
# whole cohort is 727 mothers
cohort = mn_demo %>%
left_join(., mn_outcome, by = "SID") %>%
mutate(Subsample = "Entire Cohort") %>%
rename(HOME_SCORE = HOMETOT)
demo = bind_rows(cohort, ppp_343, ppp_311, ppp_289) %>%
mutate_at(vars(SMOKER_IN_HOME, SEX, ALCOHOL), as_factor) %>%
dplyr::select(ETH, M_AGE, M_EDU, M_IQ, HOME_SCORE, mat_hard, MARITAL_STATUS, ALCOHOL, SMOKER_IN_HOME,
SEX, WISC, Subsample) %>%
rename(Ethnicity = ETH,
`Maternal education` = M_EDU,
`Marital status` = MARITAL_STATUS,
`Child sex` = SEX,
`HOME scale` = HOME_SCORE,
`Maternal IQ` = M_IQ,
`Maternal age at delivery` = M_AGE,
`Prenatal alcohol consumption` = ALCOHOL,
`Smoker in home` = SMOKER_IN_HOME,
`Material hardship` = mat_hard,
`Child Full Scale IQ` = WISC) %>%
mutate(Subsample = fct_inorder(Subsample))
Table 1
demo %>%
filter(Subsample == "BN2MF Subsample") %>%
tbl_summary(missing = "no",
statistic = list(all_continuous() ~ "{mean} ({sd})")) %>%
as_gt() %>%
tab_header(title = "Subject demographics and distribution of potential confounders, model covariates, and outcome variable")
| Characteristic |
N = 343 |
| Ethnicity |
|
| African American |
119 (35%) |
| Hispanic/Latina |
224 (65%) |
| Maternal age at delivery |
25.6 (4.9) |
| Maternal education |
|
| ≥ High school degree or equivalent |
216 (63%) |
| < High school degree |
127 (37%) |
| Maternal IQ |
85 (14) |
| HOME scale |
39.5 (6.3) |
| Material hardship |
144 (42%) |
| Marital status |
|
| Never Married |
232 (68%) |
| Ever Married |
111 (32%) |
| Prenatal alcohol consumption |
91 (27%) |
| Smoker in home |
103 (30%) |
| Child sex |
|
| Male |
184 (54%) |
| Female |
159 (46%) |
| Child Full Scale IQ |
97 (13) |
| Subsample |
|
| Entire Cohort |
0 (0%) |
| BN2MF Subsample |
343 (100%) |
| WISC Subsample |
0 (0%) |
| Complete Cases |
0 (0%) |
# demo %>% filter(Subsample == "Study Population") %>%
# tbl_summary(missing = "no", digits = list(all_continuous() ~ 1,
# all_categorical() ~ 1),
# statistic = list(all_continuous() ~ "{mean} & {sd}",
# all_categorical() ~ "{n} & {p}")) %>%
# as_kable_extra(., format = "latex")
# ^ print to latex
Supplemental table 1
demo %>%
tbl_summary(by = Subsample, missing = "no",
statistic = list(all_continuous() ~ "{mean} ({sd})")) %>% add_p()
| Characteristic |
Entire Cohort, N = 727 |
BN2MF Subsample, N = 343 |
WISC Subsample, N = 311 |
Complete Cases, N = 289 |
p-value |
| Ethnicity |
|
|
|
|
>0.9 |
| African American |
254 (35%) |
119 (35%) |
107 (34%) |
100 (35%) |
|
| Hispanic/Latina |
473 (65%) |
224 (65%) |
204 (66%) |
189 (65%) |
|
| Maternal age at delivery |
25.2 (4.9) |
25.6 (4.9) |
25.5 (4.9) |
25.6 (4.9) |
0.4 |
| Maternal education |
|
|
|
|
>0.9 |
| ≥ High school degree or equivalent |
456 (64%) |
216 (63%) |
196 (63%) |
186 (64%) |
|
| < High school degree |
257 (36%) |
127 (37%) |
115 (37%) |
103 (36%) |
|
| Maternal IQ |
85 (13) |
85 (14) |
85 (13) |
85 (13) |
0.7 |
| HOME scale |
39.4 (6.3) |
39.5 (6.3) |
39.3 (6.3) |
39.4 (6.2) |
>0.9 |
| Material hardship |
321 (44%) |
144 (42%) |
128 (41%) |
116 (40%) |
0.6 |
| Marital status |
|
|
|
|
0.9 |
| Never Married |
473 (65%) |
232 (68%) |
209 (67%) |
195 (67%) |
|
| Ever Married |
250 (35%) |
111 (32%) |
102 (33%) |
94 (33%) |
|
| Prenatal alcohol consumption |
192 (26%) |
91 (27%) |
77 (25%) |
72 (25%) |
>0.9 |
| Smoker in home |
246 (34%) |
103 (30%) |
94 (30%) |
89 (31%) |
0.5 |
| Child sex |
|
|
|
|
>0.9 |
| Male |
376 (52%) |
184 (54%) |
166 (53%) |
155 (54%) |
|
| Female |
351 (48%) |
159 (46%) |
145 (47%) |
134 (46%) |
|
| Child Full Scale IQ |
99 (13) |
97 (13) |
97 (13) |
97 (13) |
0.4 |
# demo %>%
# tbl_summary(by = Subsample,
# missing = "no", digits = list(all_continuous() ~ 1,
# all_categorical() ~ 1),
# statistic = list(all_continuous() ~ "{mean} & {sd}",
# all_categorical() ~ "{n} & {p}")) %>% add_p() %>%
# as_kable_extra(., format = "latex")
# ^^ supplemental table
Distributions & dectection
# mn_phenol and mn_pht are datasets I have stored internally in a package
lod <- mn_phenol %>% drop_na(BP_3:BPA) %>%
rename(B_PB_detect = B_PB_LOD_0_2, TCS_detect = TCS_LOD_2_3,
DCP_24_detect = X24_DCP_LOD_0_2, BP_3_detect = BP_3_LOD_0_4) %>%
mutate_at(vars(10:13), ~ifelse(is.na(.), 0, 1)) %>% # 1 if detected, 0 otherwise
inner_join(., mn_pht, by = "SID") %>%
dplyr::select(grep("detect", colnames(.))) %>%
gather(chem, level) %>%
mutate(level = ifelse(level > 1, 0, level)) %>%
group_by(chem, level) %>%
summarize(n = n()) %>%
spread(level, n) %>%
mutate_all(~replace_na(., 0))
# Number of obs with each flag for each pht
# 14 bc 3 detected 100%
detected <- lod %>%
mutate(Prop_Detected = `0`/(`1` + `0`),
`% <LOD` = 1 - Prop_Detected) %>%
ungroup() %>%
dplyr::select(Chemical = chem, `% <LOD`) %>%
mutate(Chemical = str_remove(Chemical, "_detect")) %>%
rbind(., tibble(Chemical = "M_PB", `% <LOD` = 0)) %>% # add the three detected 100%
rbind(., tibble(Chemical = "P_PB", `% <LOD` = 0)) %>%
rbind(., tibble(Chemical = "DCP_25", `% <LOD` = 0))
conc = mn_ppp %>% # get summary stats for concentrations
pivot_longer(MEHHP:BPA,
names_to = "Chemical") %>%
group_by(Chemical) %>%
summarise(Mean = mean(value, na.rm = TRUE),
Stdev = sd(value, na.rm = TRUE),
Min = min(value, na.rm = TRUE),
qs = quantile(value, c(0.25, 0.5, 0.75), na.rm = TRUE), prob = c("Q25", "Median", "Q75"),
Max = max(value, na.rm = TRUE)) %>%
pivot_wider(names_from = "prob",
values_from = "qs") %>%
dplyr::select(Chemical:Min, Q25:Q75, Max)
ppp_table <- left_join(detected, conc, by = c("Chemical")) %>% # combine concentrations and detection into table
mutate(`% <LOD` = `% <LOD` * 100) %>%
mutate_if(is.numeric, round, digits = 4) %>%
mutate(group = case_when(Chemical == "M_PB" ~ "phenols",
grepl("^M", Chemical) ~ "phthalates",
TRUE ~ "phenols")) %>%
arrange(group, Chemical) %>%
dplyr::select(-group)
Table 2
ppp_table %>% knitr::kable(caption = "Distribution of phathlate metabolites & phenols (ng/ml) in maternal spot urine during the third trimester of pregnancy (n = 343)")
Distribution of phathlate metabolites & phenols (ng/ml) in maternal spot urine during the third trimester of pregnancy (n = 343)
| B_PB |
38.4840 |
2.7956 |
7.3824 |
0.0754 |
0.1678 |
0.3692 |
1.9739 |
74.8568 |
| BP_3 |
0.5831 |
71.4030 |
215.1615 |
0.5657 |
4.7636 |
8.9000 |
26.2812 |
1999.9959 |
| BPA |
6.7055 |
2.8727 |
4.1246 |
0.2263 |
1.1038 |
1.8462 |
3.1750 |
47.1999 |
| DCP_24 |
0.5831 |
8.7942 |
17.0101 |
0.4000 |
1.6000 |
2.9714 |
6.0439 |
114.2852 |
| DCP_25 |
0.0000 |
206.4596 |
306.9856 |
4.4445 |
37.9913 |
81.8461 |
205.5894 |
2003.2019 |
| M_PB |
0.0000 |
256.3146 |
318.7391 |
3.2000 |
46.0640 |
132.4444 |
363.6577 |
2026.6525 |
| P_PB |
0.0000 |
72.2174 |
128.6541 |
0.4571 |
5.3333 |
20.4000 |
71.9046 |
941.1787 |
| TCS |
19.8251 |
71.2825 |
164.0853 |
0.8973 |
3.5200 |
8.6738 |
36.6539 |
1080.0032 |
| MBP |
0.0000 |
55.5034 |
59.7337 |
4.0800 |
22.7634 |
37.0824 |
66.4067 |
657.7790 |
| MBZP |
0.0000 |
28.3870 |
49.6796 |
0.8436 |
6.5874 |
13.0226 |
27.3599 |
419.3823 |
| MCPP |
4.3732 |
2.9842 |
2.7078 |
0.0943 |
1.4000 |
2.3429 |
3.7112 |
26.6999 |
| MECPP |
0.0000 |
73.3590 |
138.7470 |
4.5091 |
21.0476 |
35.3883 |
69.0421 |
1382.8514 |
| MEHHP |
0.0000 |
51.5927 |
124.3645 |
0.8727 |
11.2500 |
20.9778 |
40.0445 |
1217.3912 |
| MEHP |
15.7434 |
13.2314 |
33.5140 |
0.5903 |
2.2797 |
4.9185 |
11.5022 |
426.4347 |
| MEOHP |
0.0000 |
39.3903 |
88.8139 |
1.0182 |
9.7615 |
17.0666 |
32.5264 |
918.2608 |
| MEP |
0.0000 |
335.3467 |
597.1151 |
14.0801 |
77.5560 |
143.6121 |
333.7031 |
7440.7320 |
| MIBP |
0.5831 |
13.6081 |
18.5089 |
0.6857 |
5.4978 |
9.6000 |
16.0785 |
272.2915 |
# ppp_table %>%
# dplyr::select(1:2, 6:8) %>%
# stargazer::stargazer(summary = F)
# ^ print to latex
ewa_raw %>%
summarise_all(quantile, probs = seq(0.25, 0.75, 0.25)) %>% bind_cols(probs = seq(0.25, 0.75, 0.25), .)
| probs | P1 | P2 |
| 0.25 | 2.85 | 2.18 |
| 0.5 | 4.21 | 3.41 |
| 0.75 | 6.13 | 6.91 |
# added this to table
Chemicals
mn_ppp %>%
pivot_longer(MEHHP:BPA) %>%
mutate(name = str_replace(name, "_", "-")) %>%
ggplot(aes(x = value)) +
geom_histogram(aes(y = ..density..), alpha=0.65, fill = "#0072B5") +
facet_wrap(name~., scales = "free") +
labs(y = "Density", x="") +
theme(axis.text.y = element_blank(),
axis.text.x = element_text(size = 10))

Correlation
cormat <- get_lower_tri(round(cor(mn_ppp[,-1], use = "complete.obs",
method = c("spearman")),2))[2:17, 1:16]
melted_cormat <- melt(cormat) %>%
rename(Correlation = value) %>%
mutate(Var1 = case_when(Var1 == "DCP_24" ~ "24-DCP",
Var1 == "DCP_25" ~ "25-DCP",
Var1 == "BP_3" ~ "BP-3",
TRUE ~ as.character(Var1)),
Var2 = case_when(Var2 == "DCP_24" ~ "24-DCP",
Var2 == "DCP_25" ~ "25-DCP",
Var2 == "BP_3" ~ "BP-3",
TRUE ~ as.character(Var2)),
Var1 = fct_inorder(str_remove(Var1, "_")),
Var2 = fct_inorder(str_remove(Var2, "_"))) %>% as_tibble()
Exposure patterns
Loadings
applied <- eh %>%
as_tibble() %>%
mutate(Pattern = 1:2) %>%
gather(key = Chemicals, value = Loadings, -Pattern) %>%
mutate(Class = as_factor(case_when(grepl("PB", Chemicals) ~ "Parabens",
grepl("^M[A-Z]", Chemicals) ~ "Phthalates",
TRUE ~ "Phenols")),
Chemicals = str_replace(Chemicals, "_", "-")) %>%
mutate(Pattern = paste0("Pattern ", Pattern),
Chemicals = fct_inorder(Chemicals))
Health models
Unadjusted
fit_p1 <- lm(WISC ~ P1 + SEX + SEX*P1, data = mn)
add_ci4interaction(fit_p1, "P1", "P1:SEXFemale")
| term | estimate | std.error | statistic | p.value | 2.5 % | 97.5 % |
| (Intercept) | 98 | 1.66 | 58.9 | 1.04e-168 | 94.8 | 101 |
| P1 in males | -0.0673 | 1.08 | -0.0626 | 0.95 | -2.18 | 2.05 |
| SEXFemale | 2.1 | 2.69 | 0.778 | 0.437 | -3.21 | 7.4 |
| P1:SEXFemale | -3.32 | 1.94 | -1.72 | 0.0873 | -7.14 | 0.489 |
| P1 in females | -3.39 | 1.61 | 2.1 | 0.0351 | -6.55 | -0.232 |
fit_p2 <- lm(WISC ~ P2 + SEX + SEX*P2, data = mn)
add_ci4interaction(fit_p2, "P2", "P2:SEXFemale")
| term | estimate | std.error | statistic | p.value | 2.5 % | 97.5 % |
| (Intercept) | 100 | 1.75 | 57.3 | 2.07e-165 | 96.7 | 104 |
| P2 in males | -1.55 | 1 | -1.55 | 0.123 | -3.52 | 0.422 |
| SEXFemale | -6.87 | 2.59 | -2.65 | 0.00847 | -12 | -1.77 |
| P2:SEXFemale | 3.86 | 1.55 | 2.48 | 0.0135 | 0.801 | 6.91 |
| P2 in females | 2.31 | 1.19 | 1.94 | 0.0514 | -0.0181 | 4.63 |
Adjusted
Pattern 1
fit_p1a <- lm(WISC ~ P1 + SEX*P1 + SEX + M_IQ + ALCOHOL + M_EDU + M_AGE +
MARITAL_STATUS + HOME_SCORE + mat_hard
, data = mn)
add_ci4interaction(fit_p1a, "P1", "P1:SEXFemale")
| term | estimate | std.error | statistic | p.value | 2.5 % | 97.5 % |
| (Intercept) | 102 | 2.06 | 49.4 | 4.47e-139 | 97.6 | 106 |
| P1 in males | -0.363 | 1.03 | -0.351 | 0.726 | -2.4 | 1.67 |
| SEXFemale | 1.3 | 2.62 | 0.497 | 0.62 | -3.85 | 6.46 |
| M_IQ | 3.12 | 0.813 | 3.84 | 0.000155 | 1.52 | 4.72 |
| ALCOHOLYes | -5.45 | 1.7 | -3.2 | 0.00151 | -8.8 | -2.1 |
| M_EDU< High school degree | -1.29 | 1.65 | -0.782 | 0.435 | -4.55 | 1.96 |
| M_AGE | 1.06 | 0.845 | 1.25 | 0.212 | -0.605 | 2.72 |
| MARITAL_STATUSEver Married | -0.146 | 1.71 | -0.0855 | 0.932 | -3.5 | 3.21 |
| HOME_SCORE | 2.34 | 0.775 | 3.02 | 0.00277 | 0.814 | 3.86 |
| mat_hardYes | -2.62 | 1.53 | -1.71 | 0.0887 | -5.63 | 0.399 |
| P1:SEXFemale | -3.08 | 1.86 | -1.66 | 0.0987 | -6.75 | 0.581 |
| P1 in females | -3.45 | 1.55 | 2.22 | 0.0265 | -6.49 | -0.398 |
glance(fit_p1a)[,-c(11:12)]
| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC | BIC | deviance |
| 0.177 | 0.147 | 12.4 | 5.93 | 4e-08 | 10 | -1.12e+03 | 2.27e+03 | 2.31e+03 | 4.22e+04 |
Pattern 2
fit_p2a <- lm(WISC ~ P2 + SEX*P2 + SEX + M_IQ + ALCOHOL + M_EDU + M_AGE +
MARITAL_STATUS + HOME_SCORE + mat_hard
, data = mn)
add_ci4interaction(fit_p2a, "P2", "P2:SEXFemale")
| term | estimate | std.error | statistic | p.value | 2.5 % | 97.5 % |
| (Intercept) | 102 | 2 | 51 | 1.63e-142 | 98.3 | 106 |
| P2 in males | -0.803 | 0.982 | -0.817 | 0.414 | -2.74 | 1.13 |
| SEXFemale | -5.99 | 2.56 | -2.34 | 0.02 | -11 | -0.952 |
| M_IQ | 2.93 | 0.81 | 3.62 | 0.000347 | 1.34 | 4.53 |
| ALCOHOLYes | -5.54 | 1.71 | -3.24 | 0.00133 | -8.91 | -2.18 |
| M_EDU< High school degree | -1.27 | 1.66 | -0.766 | 0.444 | -4.53 | 1.99 |
| M_AGE | 0.96 | 0.85 | 1.13 | 0.26 | -0.713 | 2.63 |
| MARITAL_STATUSEver Married | -0.153 | 1.71 | -0.0898 | 0.929 | -3.52 | 3.21 |
| HOME_SCORE | 2.37 | 0.777 | 3.05 | 0.0025 | 0.842 | 3.9 |
| mat_hardYes | -2.27 | 1.52 | -1.49 | 0.137 | -5.27 | 0.728 |
| P2:SEXFemale | 2.8 | 1.51 | 1.85 | 0.0656 | -0.182 | 5.78 |
| P2 in females | 2 | 1.15 | 1.74 | 0.0814 | -0.252 | 4.25 |
glance(fit_p2a)[,-c(11:12)]
| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC | BIC | deviance |
| 0.173 | 0.143 | 12.4 | 5.77 | 7.11e-08 | 10 | -1.12e+03 | 2.27e+03 | 2.32e+03 | 4.24e+04 |
Table 3
tabp1 = add_ci4interaction(fit_p1a, "P1", "P1:SEXFemale")[c(2,11:12),] %>% mutate(pattern = "Pattern 1")
tabp2 = add_ci4interaction(fit_p2a, "P2", "P2:SEXFemale")[c(2,11:12),] %>% mutate(pattern = "Pattern 2")
model_table = bind_rows(tabp1, tabp2)[,c(1:2, 6:8)] %>% mutate(model = "Traditional") %>%
dplyr::select(model, pattern, everything()) %>%
arrange(model, pattern, term) %>%
mutate_if(is.numeric, round, 1) %>%
mutate(`95% Confidence Interval` = str_c("(", `2.5 %`, ", ", `97.5 %`, ")")) %>%
mutate(term = case_when(term == "sex*p" | grepl("SEX", term) ~ " interaction term",
TRUE ~ term),
term = str_remove(term, "(P1|P2|p)"),
term = str_c(pattern, term)) %>%
dplyr::select(-`2.5 %`, -`97.5 %`, -pattern)
model_table
| model | term | estimate | 95% Confidence Interval |
| Traditional | Pattern 1 in females | -3.4 | (-6.5, -0.4) |
| Traditional | Pattern 1 in males | -0.4 | (-2.4, 1.7) |
| Traditional | Pattern 1 interaction term | -3.1 | (-6.7, 0.6) |
| Traditional | Pattern 2 in females | 2 | (-0.3, 4.2) |
| Traditional | Pattern 2 in males | -0.8 | (-2.7, 1.1) |
| Traditional | Pattern 2 interaction term | 2.8 | (-0.2, 5.8) |
# health model table
# stargazer::stargazer(model_table, summary = F)